Redlining Analysis: CES Pollution Burden - Environmental Effects
Updated: 2020-05-05
# Markdown
# library(bookdown)
# General data analysis and transformation
library(readr)
library(readxl)
library(dplyr)
library(janitor)
# library(gtools)
library(stringr)
library(units)
library(tidyr)
library(rlang)
# API-related
# library(jsonlite)
# library(urltools)
# Mapping and GIS operations
library(sf)
# library(leaflet)
# library(htmlwidgets)
# library(geojsonsf)
# library(rmapshaper)
library(tmap)
library(tmaptools)
library(ceramic)
# workflow
library(here)
# plotting
library(ggplot2)
library(forcats)# Get list of CES 3 variable names and plain text names
ces3_variables <- readr::read_csv(here('data_prepared', 'ces_names.csv'))
ces_measures_effects <- ces3_variables %>%
filter(group == 'pollution burden',
subgroup == 'environmental effects',
type == 'percentile')
ces_measures_exposures <- ces3_variables %>%
filter(group == 'pollution burden',
subgroup == 'exposures',
type == 'percentile')
ces_measures_sensitive_pops <- ces3_variables %>%
filter(group == 'population characteristics',
subgroup == 'sensitive populations',
type == 'percentile')
ces_measures_socioeconomic <- ces3_variables %>%
filter(group == 'population characteristics',
subgroup == 'socioeconomic factors',
type == 'percentile')Introduction
This document attempts to do some quantitative geospatial analysis to look at correlations between redline mapping and indicators of water quality, public health, and facility information. It computes the area weighed average percentile for different CalEnviroScreen 3.0 (CES) indicators by HOLC rating for a given city and statewide. CES considers two broad categories of indicators, (1) Pollution Burden and (2) Population Characteristics, which are both further broken down into two sub-categories each (factsheet, indicators page).
The Pollution Burden category includes the following indicators:
- Environmental Effects
- Cleanups Percentile
- GroundWater Threats Percentile
- Haz Waste Percentile
- Imp Water Bodies Percentile
- Solid Waste Percentile
- Exposures
- Ozone Percentile
- PM 2.5 Percentile
- Diesel PM Percentile
- Drinking Water Percentile
- Pesticides Percentile
- Tox Releases Percentile
- Traffic Percentile
The Population Characteristics category includes the following indicators:
- Sensitive Populations
- Asthma Percentile
- Low Birth Weight Percentile
- Cardiovascular Disease Percentile
- Socioeconomic Factors
- Education Percentile
- Linguistic Isolation Percentile
- Poverty Percentile
- Unemployment Percentile
- Housing Burden Percentile
Analysis Notes
The CES polygons and the Redline polygons have different coverages, as shown in the maps below. In the first map for each indicator, the left pane of each map shows the HOLC rated areas, the center pane shows the CES indicator scores (at the census tract level) with the Redlined areas superimposed on top, and the right pane shows the overlapping portions of the CES polygons and the redlined areas (border colors represent HOLC rating).
# load data
# Redline Polygons
redline_polygons <- read_rds(here('data_prepared', 'redline_polygons.RDS'))
# st_crs(redline_polygons) # to check the reference system
# modify the 'city' column to better keep track when joining with other datasets
redline_polygons <- redline_polygons %>% rename('redline_city' = 'city')
# add the area of each redline polygon
redline_polygons <- redline_polygons %>% mutate(redline_poly_area = st_area(.))
# NOTE: not all redline polygons have an HOLC ID listed - create a unique ID for each one with a missing ID (so that each unique HOLC polygon can be associated with a weighted average CES score)
# create new column for the new holc id to be assigned
redline_polygons <- redline_polygons %>%
mutate('holc_id_2' = holc_id) %>%
arrange(redline_city, holc_grade, holc_id)
# initialiaze some variables for the loop below
grade_city_vec <- c()
j <- 1
# loop through all HOLC polygons ...
for (i_row in 1:nrow(redline_polygons)) {
cur_id <- redline_polygons$holc_id[i_row]
grade <- redline_polygons$holc_grade[i_row]
city <- redline_polygons$redline_city[i_row]
# the current combination of grade and city
grade_city_cur <- paste0(city, '_', grade)
# if the holc id is NA, check to see if there are already other missing holc id's for that city/grade combination - if so, increment a counter by 1, otherwise set counter (back) to 1
if (is.na(cur_id)) {
if(grade_city_cur %in% grade_city_vec) {
grade_city_counter <- grade_city_counter + 1
} else {
grade_city_counter <- 1
}
# enter a new unique ID for the row with the missing ID
redline_polygons$holc_id_2[i_row] <- paste0('unknown_',
grade, '_',
grade_city_counter)
# add the grade/city combo to a running list to track number of missing holc id's for that city/grade combination
grade_city_vec[j] <- grade_city_cur
j <- j + 1
}
}
# CES 3 Polygons
ces_3_poly <- read_sf(here('data_raw', 'CalEnviroScreen3', 'CESJune2018Update_SHP','CES3June2018Update.shp'))
# st_crs(ces_3_poly) # to check the reference system
# modify the 'City' column to better keep track when joining with other datasets
ces_3_poly <- ces_3_poly %>% rename('CES_City' = 'City')
# Transform both to UTM Zone 10N
redline_polygons <- redline_polygons %>% st_transform(crs = 26910)
ces_3_poly <- ces_3_poly %>% st_transform(crs = 26910)
# create a map to check that the data loaded correctly
# tm_shape(ces_3_poly) + tm_borders() + tm_shape(redline_polygons) + tm_borders(col = 'red') # get a basemap (See: https://github.com/mtennekes/tmap/issues/185)
Sys.setenv(MAPBOX_API_KEY = "pk.eyJ1IjoiZGFsdGFyZSIsImEiOiJjazZqcm5ocHgwMG84M2tteGc5eWd4YXNmIn0.nHiYenoq-wFWwH3V4vtPXA")
background <- cc_location(ces_3_poly %>% filter(CES_City == 'Sacramento'), verbose = FALSE) # create map
map_redline <- tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
tm_polygons('holc_grade', palette = c('green', 'blue', 'yellow', 'red'), title = 'HOLC Grade') +
tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE) # get the intersection of the CES polygons and redline polygons
ces3_redline_intersection <- st_intersection(ces_3_poly, redline_polygons)
# inspect
# glimpse(ces3_redline_intersection)
# add the area of each clipped polygon to the data frame
ces3_redline_intersection <- ces3_redline_intersection %>%
mutate(clipped_area = st_area(.))
# glimpse(ces3_redline_intersection)
# # write to output file
# st_write(obj = ces3_redline_intersection,
# here('analysis-output', 'ces3_redline_intersection.shp')) ces_poly_map <- function(ces_id, ces_measure_name) {
# plot the CES polygons
map_ces <- tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento' | CES_City == 'West Sacramento')) +
tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 1) +
tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
tm_borders(lwd = 0, alpha = 0) +
# tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
# tm_borders(lwd = 2, col = 'black') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'D')) +
tm_borders(lwd = 2, col = 'red') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'C')) +
tm_borders(lwd = 2, col = 'gold2') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'B')) +
tm_borders(lwd = 2, col = 'blue') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'A')) +
tm_borders(lwd = 2, col = 'green') +
# tm_text('holc_grade', size = 0.5) +
tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
map_ces
} ces_redline_overlap_map <- function(ces_id, ces_measure_name) {
map_overlap <- tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento')) +
tm_borders(lwd = 1) +
tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) +
# tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
# tm_borders(lwd = 2, col = 'black') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'D')) +
tm_borders(lwd = 2, col = 'red') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'C')) +
tm_borders(lwd = 2, col = 'gold2') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'B')) +
tm_borders(lwd = 2, col = 'blue') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'A')) +
tm_borders(lwd = 2, col = 'green') +
#tm_text('holc_grade', size = 0.5) +
tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
map_overlap
} # FACET MAP
ces_redline_facet_map <- function(ces_id, ces_measure_name) {
map_facets <- tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento')) +
tm_borders(lwd = 1) +
tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(ces3_redline_intersection %>%
filter(redline_city == 'Sacramento'), is.master = TRUE) +
# tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento')) +
tm_facets(by = 'holc_grade', free.coords = FALSE) +
tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) +
tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
tm_facets(by = 'holc_grade') +
tm_borders(lwd = 2, col = 'black') +
# # D
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'D')) +
# tm_facets(by = 'holc_grade') +
# tm_borders(lwd = 2, col = 'red') +
# # C
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'C')) +
# tm_borders(lwd = 2, col = 'gold2') +
# tm_facets(by = 'holc_grade') +
# # B
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'B')) +
# tm_borders(lwd = 2, col = 'blue') +
# tm_facets(by = 'holc_grade') +
# # A
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'A')) +
# tm_borders(lwd = 2, col = 'green') +
# tm_facets(by = 'holc_grade') +
#
tm_legend(bg.color= 'white', frame = TRUE)
map_facets
}ces_scores_plot <- function(ces_id, ces_measure_name) {
# group by city and holc rating, then compute the weighted average score for each
ces3_redline_grouped <- ces3_redline_intersection %>%
st_drop_geometry() %>%
mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>%
group_by(redline_city, holc_grade) %>%
summarize(total_area = sum(clipped_area),
area_x_score_total = sum(area_x_score)) %>%
mutate(weighted_score = area_x_score_total / total_area) %>%
drop_units()
ces3_redline_grouped_statewide <- ces3_redline_grouped %>%
group_by(holc_grade) %>%
summarize(total_area = sum(total_area),
area_x_score_total = sum(area_x_score_total)) %>%
mutate(weighted_score = area_x_score_total / total_area) %>%
mutate(redline_city = '**Statewide**') %>%
drop_units()
ces3_redline_grouped_city_state <- bind_rows(ces3_redline_grouped, ces3_redline_grouped_statewide)
# plot the results by city and holc rating, and include statewide summary
g_city_state <- ggplot(ces3_redline_grouped_city_state) +
aes(x = fct_relevel(fct_reorder(redline_city, weighted_score),
c('**Statewide**'),
after = 0),
y = weighted_score,
fill = fct_rev(holc_grade)) +
geom_bar(stat = 'identity', position = 'dodge') +
geom_vline(xintercept = 1.52, size = 0.5, color = 'grey50') +
scale_fill_manual(values = rev(c('green', 'blue', 'yellow', 'red'))) +
guides(fill = guide_legend(reverse = TRUE)) +
labs(fill = 'HOLC Rating',
title = ces_measure_name,
x = NULL,
y = "Area Weighted Average Score By City") +
coord_flip()
g_city_state
}Pollution Burden - Environmental Effects
Cleanups Percentile
i_map <- 1
ces_measure_name <- ces_measures_effects$name[i_map]
ces_id <- ces_measures_effects$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
GroundWater Threats Percentile
i_map <- 2
ces_measure_name <- ces_measures_effects$name[i_map]
ces_id <- ces_measures_effects$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Haz Waste Percentile
i_map <- 3
ces_measure_name <- ces_measures_effects$name[i_map]
ces_id <- ces_measures_effects$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Imp Water Bodies Percentile
i_map <- 4
ces_measure_name <- ces_measures_effects$name[i_map]
ces_id <- ces_measures_effects$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Solid Waste Percentile
i_map <- 5
ces_measure_name <- ces_measures_effects$name[i_map]
ces_id <- ces_measures_effects$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Pollution Burden - Exposures
Ozone Percentile
i_map <- 1
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
PM 2.5 Percentile
i_map <- 2
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Diesel PM Percentile
i_map <- 3
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Drinking Water Percentile
i_map <- 4
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Pesticides Percentile
i_map <- 5
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Tox Releases Percentile
i_map <- 6
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Traffic Percentile
i_map <- 7
ces_measure_name <- ces_measures_exposures$name[i_map]
ces_id <- ces_measures_exposures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Population Characteristics - Sensitive Populations
Asthma Percentile
i_map <- 1
ces_measure_name <- ces_measures_sensitive_pops$name[i_map]
ces_id <- ces_measures_sensitive_pops$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Low Birth Weight Percentile
i_map <- 2
ces_measure_name <- ces_measures_sensitive_pops$name[i_map]
ces_id <- ces_measures_sensitive_pops$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Cardiovascular Disease Percentile
i_map <- 3
ces_measure_name <- ces_measures_sensitive_pops$name[i_map]
ces_id <- ces_measures_sensitive_pops$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Population Characteristics - Socioeconomic Factors
Education Percentile
i_map <- 1
ces_measure_name <- ces_measures_socioeconomic$name[i_map]
ces_id <- ces_measures_socioeconomic$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Linguistic Isolation Percentile
i_map <- 2
ces_measure_name <- ces_measures_socioeconomic$name[i_map]
ces_id <- ces_measures_socioeconomic$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Poverty Percentile
i_map <- 3
ces_measure_name <- ces_measures_socioeconomic$name[i_map]
ces_id <- ces_measures_socioeconomic$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Unemployment Percentile
i_map <- 4
ces_measure_name <- ces_measures_socioeconomic$name[i_map]
ces_id <- ces_measures_socioeconomic$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Housing Burden Percentile
i_map <- 5
ces_measure_name <- ces_measures_socioeconomic$name[i_map]
ces_id <- ces_measures_socioeconomic$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
Create Output Shapefile
The code below creates a shapefile with the calculated weighted average scores for all HOLC polygons across all of the CES indicators.
Note that the field names are the same as those used in the original CES shapefile, but in the output shapefile produced by the code below the values refer to the calculated weighted average scores.
Also, note that this output shapefile contains a field called ‘holc_id_2’ rather than ‘holc_id’. This ‘holc_id_2’ is a revised version of the ‘holc_id_field’ - this is needed because not all HOLC polygons are labeled with an HOLC ID in the source shapefiles (from the University of Richmond site), so it was necessary to create a unique HOLC ID within each city/grade combination so that the clipped CES polygons could all be grouped to a unique HOLC polygon. In the ‘holc_id_2’ field, those records with missing HOLC IDs have an entry with the format “unknown_X_Y”, where X is the HOLC grade and Y is a number.
ces_measures_combined <- bind_rows(ces_measures_effects,
ces_measures_exposures,
ces_measures_sensitive_pops,
ces_measures_socioeconomic)
for (i_output in seq(nrow(ces_measures_combined))) {
ces_measure_name <- ces_measures_combined$name[i_output]
ces_id <- ces_measures_combined$id[i_output]
var_name <- ces_id # paste0('weighted_score_', ces_id)
# group by city and holc ID, then compute the weighted average score for each holc ID
ces3_redline_grouped_holc_id <- ces3_redline_intersection %>%
mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>%
group_by(redline_city, holc_id_2, holc_grade) %>%
summarize(total_area = sum(clipped_area),
area_x_score_total = sum(area_x_score),
component_poly_count = n()) %>%
mutate(!!var_name := area_x_score_total / total_area) %>%
drop_units()
ces3_redline_grouped_holc_id <- ces3_redline_grouped_holc_id %>%
select(redline_city, holc_id_2, holc_grade, !!var_name) %>%
rename(city = redline_city)
# add to output df
if (i_output == 1) {
df_output <- ces3_redline_grouped_holc_id
} else {
df_output <- left_join(x = df_output,
y = ces3_redline_grouped_holc_id %>%
st_drop_geometry(),
by = c('city', 'holc_id_2', 'holc_grade'))
}
}
df_output <- df_output %>% arrange(city, holc_grade, holc_id_2)
# write the output shapefile
st_write(obj = df_output,
here('analysis-output',
'redline_ces_scores',
paste0('redline_ces_scores',
'_',
Sys.Date(),
'.shp'))
)## Writing layer `redline_ces_scores_2020-05-05' to data source `C:/David/Redline-Mapping/analysis-output/redline_ces_scores/redline_ces_scores_2020-05-05.shp' using driver `ESRI Shapefile'
## Writing 868 features with 23 fields and geometry type Unknown (any).